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ABSTRACT 


A numerical method for solving a problem in unsteady 
slag flow in the hearth of a blast furnace is presented. 
This problem is reduced to some free boundary problem 
for an elliptic system. The potential problem for a 
given free boundary is approximated by the penalty meth¬ 
od. The derivatives of the potential function on the 
free boundary is approximated by the integration of the 
penalty term, and then the subsequent shape of the free 
boundary is obtained by solving the differential equa¬ 
tion for the motion of the free boundary. The finite 
difference method is used to solve the penalized prob¬ 
lem. A numerical example is given. 
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Numerical Solution of Free Boundary Problem 
for Unsteady Slag Flow in the Hearth 


Makoto NATORI and Hideo KAWARADA 


1. Introduction 

We present a numerical method for solving a problem in unsteady 
flow of molton slag in the hearth region of iron producing blast 
furnaces during the tapping operation[1]. This problem is reduced 
to v some free boundary problems for an elliptic system. This type 
of problem is similar to the porous flow of underground water in 
which the water surface is a free boundary. The numerical calcula¬ 
tions of this type of problem were done by various researchers[2-4]. 
The three-dimensional problem of the slag flow in the hearth was 
solved by using the finite element method by Ichihara and Fukutake[5]. 
They concluded that their computation scheme is not so efficient as 
it is applied for the practical use. 

The object of this paper is to settle this computational insta¬ 
bility by using the penalty method developed by Kawarada and Natori 
[6-9]. 


1 



2. Formulation 


We consider two-dimensional slag flow in the hearth which is 
bounded by impermeable boundaries y=0, x=0 and x=a. One of vertical 
boundaries, x=0, has a tapping hole near the bottom. While y=g(x,t) 
denotes the free surface of the slag region(See Fig.l). 



A velocity potential can be defined by 
P 

4> = - + y 

Y 

where p is the fluid pressure and Y is the specific weight of the 
fluid. If it is assumed that Darcy's law holds, the potential is 
given by 

(1) A <f> = 0 

(2) 4) = y on y = g(x,t) 
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(3) <f> = 0 on y = 0 

(4) <j> x = 0 on x = 0 and x = a, except on the tapping hole 

(5) <J> = k ( >0 ) on the tapping hole. 


where k is some constant. 

The motion of the free surface is given by 


(6) g t = 


( <J) - d> g 

T y Y x ^x 




3n 


y=g(x,t) 
y=g (x,t) 


The initial shape of the free surface, y = g(x,0), is given and 
forms an initial condition for equation (6). 

When we try to solve the problem formulated above, a numerical 
procedure must contain a routine for solving the potential problem 
(1) - (5) for a given free boundary y = g(x,t). When this is done, 
the derivatives of the potential function can be calculated on the 
free boundary, and then the equation (6) can be solved for the sub¬ 
sequent shape of the free boundary. 

If we use the method of the integrated penalty to solve the 
potential problem, then the derivatives of the potential function 
on the free boundary are easily obtained. This is the reason of 
our application of the penalty method to the free boundary problems. 


3. Penalty Method 
3.1 Penalized problem 

We define the characteristic function x E (x,y,t) such as 


(7) x £ (x,y,t) 


1 in fl-ft e 

g 

0 in ft e 

g 

— 3 — 



where the domains ft e and ft, which includes fie, are defined by 
ft e = { (x ,y) 10 < x < a, 0 < y < g £ (x,t)} 

and 

ft = {(x,y)|0 <x<a, 0<y<b} 

as shown in Fig.2. 



Fig. 2 

By the use of y £ r equation (1) is approximated by 

(8) A<j> £ - |x e (* E -y) =0 in S2, 

where £ is a positive constant. We add a new boundary condition: 

(9) <\> e = y on y = b, 

to the boundary conditions (3)-(5). 

In fact, if we let e be sufficiently small then we know that 
<f) £ approximates <j> in ft e and (J) £ is nearly equal to y in ft-ft e. 
Therefore the boundary condition (2) is approximately satisfied. 
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If we use the method of integrated penalty, equation (6) is 
approximated by 

(10) = -(1- i [ X £ (<j> £ - Y) dy ) . 

o 

Now, the second term in the right hand side of (10) is called the 
integrated penalty. 


3.2 Discretization of the penalized problem 

The penalized problem (8) with the boundary conditions (3)-(5) 
and (9) is discretized by the finite difference. Also the free 
boundary equation (10) is solved by Euler's method. The intervals 
0 £ x £ a and 0 y £ b are divided into N and M equal subinter¬ 
vals of width h. The mesh size of time is denoted by At. We use 
the following notations in the discretized system: 


x. = ih 
i 


Yj = 


t, = kAt 
k 


0 < i < N 
0 £ j £ M 
0 < k 


hjk = ♦ (x i'yj'V 

g ik =g £ (x. ( t k ) 

Xijk = x'^i'Yj'V 

We define the discretized characteristic function by 


(11 


x ijk=\ 


1 - pjk 

i + Jl 

4e ^ik 


j > [g ik /h] 
j = [g ik /h] 


j < 


5 



where [ ] denotes the Gauss symbol and 


p ik = g ik /h _ [g ik /h] - 

Then the potential function <j>. satisfies the following equations 
for any k : 

h 2 

(4+— Xj_j k ) ^ijk“^i-i j k _( ^i+i j k~^i j-1 k ^i j+1 k 

h 2 

(12) e ^ijk 

(0<i<N , 0<j£M ) 

This system of linear equations is solved by the incomplete Cholesky 
decomposition combined with conjugate gradient method [11]. 

The free boundary g^ k is obtained by 


(13) 

^i k+1 

= g ik +AtF(g 

(14) 

F<g ik> 

M 

= -i + |Z X 
j=0 


- vy ik ' - efj^ A ijk y ijk 

It should be noted that x.j_j k and are determined by g ifc . 


3.3 How to choose the penalty parameter e. 

We assume e is expressed by 

(15) e = h° 

and try to find an optimal value of a to minimize the difference of 
the right sides of equations (6) and (14). For this purpose, we 
consider a simple test problem (See Fig.3) : 

Au = 0 in 
u = 1-y on x = 0 

(16) 

u = 1-x on y = 0 

u = 0 on y = 1-x 

— 6 — 




Fig. 3 


This problem has an exact solution u = 1-x-y. Therfore, 

(17) - . / 1+g 2 , .= 2, 

V x 3n'y=g(x) ' 

where g(x) = 1-x. 

Here we construct the discretized problem (p£) of the penalized 
equation : 


(18) Au e - = 0 in 


We investigate the difference between (17) and the integrated 
penalty : 


(IP) 


M 

= ^Z 


E jV 0 Xi Fij 


by varying the value of a in (15). We get that the optimal value of 
a is 3~4 for 1/8 <h< 1/16. 
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3.4 Stability condition 

Here we study the stability condition of (13). It is well 
known that the stability condition is 

(19) At|f§|<2, 
where 

f( g) = -i+| ^ x e (g)($ £ (g)-y)dy* 

If we use the property [12] 

| (|) e (x,g (x,t) ) | £ C Q yfT, 

£ £ 

where Cq depends on (J) and g , then we have 

i — i S i —i 

| F (g) -F (g) | £ — | g-g | . 

F 

If we substitute C Q /fe into |3F/9g| in (19), then we have 

C 0 

0 < At—^ < 2. 

Therefore we may choose 

(20) At = 

4. Numerical Example 

In this section we present results for the problem (l)-(6), 
obtained by the method of integrated penalty. Data of the problem 
are as follows : 

a = 1 

b = 0 .3125 
k = 0 .625 

The tapping hole is located at (0,1/16). The initial surface is 
given by 
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g(x,0) = 0.25. 


The parameters used for the numerical calculations 
h = 1/16 
e = h 3 = 1/4096 
At = y/e = 1/64. 

In Fig.4 we show the results. 



are 
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